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Abstract 

The study of Cyber-Physical Systems requires practical tools for approximation and comparison of system 
behaviors. Existing approaches to these problems impose undue restrictions on the system's continuous and discrete 
dynamics. Metrization and simulation of controlled hybrid systems is considered here in a unified framework by 
constructing a state space metric. The metric is applied to develop a numerical simulation algorithm that converges 
uniformly to orbitally stable executions of controlled hybrid systems, up to and including Zeno events. Benchmark 
Fy | examples of novel hybrid phenomena illustrate the utility of the proposed tools. 

OO 

I. Introduction 

The ubiquity of Cyber-Physical Systems (CPS) demands tools for optimization and verification of com- 
Q munication and control architectures. Development of these tools in turn necessitates a formal framework 
for simulation and comparison of CPS behaviors. For continuous-state dynamical systems and finite-state 
automata there separately exist rich sets of tools applicable to these problems. The interaction of discrete 
transitions with continuous dynamics that is characteristic of CPS requires a new approach applicable to 
controlled hybrid systems. To address this problem, in this paper, we construct a distance metric over 
the state space of a controlled hybrid system and apply this metric to develop a provably-convergent 
numerical simulation algorithm. Our framework imposes mild restrictions on the system, enabling formal 
CN investigation of a wide range of CPS: the dynamics may be nonlinear, the continuous dynamics may be 
controlled, and multiple discrete transitions may occur simultaneously, so long as executions are orbitally 
T^j- stable. 

^ Efforts to topologize and subsequently metrize controlled hybrid systems have been significant but 
O fragmented. Nerode and Kohn [1] consider state-space topologies induced by the finite-state automaton 
^ associated with the hybrid system. We propose a metric topology over the state space of controlled 
hybrid systems, effectively metrizing the quotient topology proposed by Simic et al. fl2), as well as the 
topology constructed by Ames and Sastry [3] over the regularization proposed by Johansson et al. flU. In 
^ contrast, Tavernini [5] directly metrized the space of executions of hybrid systems, and Gokhman [6J later 
c3 demonstrated the equivalence of the resulting topology with that generated by the Skorokhod trajectory 
metric (see Chapter 6 in 0). We highlight the technical and practical limitations imposed by metrizing 



the trajectory space rather than state space in Section V-C 

The notion of regularization or relaxation of a hybrid system should not be confused with the "relax- 
ation" of hybrid inclusions described by Cai et al. (H. Since interpreting our controlled hybrid systems as 
hybrid inclusions yields singleton-valued "flow" and "jump" maps, relaxation in this sense does not yield 
a distinct hybrid system. Sanfelice and Teel [9J subsequently prove existence of approximating executions 
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for a given "simulation" of a hybrid inclusion. In this paper we consider the opposite problem of proving 
convergence of approximating simulations for a given execution of a controlled hybrid system. 

The literature on numerical simulation of deterministic hybrid systems may be broadly partitioned into 
two groups: practical algorithm development focused on obtaining high-precision estimates of discrete 
event times, and theoretical proofs of convergence for simulations of certain classes of hybrid systems. 
Practical algorithms aim to place time-steps close to discrete event times using root-finding ifTOll . |fTTfl. 021. 
Theoretical proofs of convergence have generally required restrictive assumptions. Esposito et al. 0"3l , 
for instance, apply feedback linearization to asymptotically guarantee event detection for semi-algebraic 
guards, while Paoli and Schatzman [fT4l develop a provably-convergent simulation algorithm for second- 
order mechanical states undergoing impact specified by a unilateral constraint. The most general con- 
vergence results relax the requirement that discrete transition times be determined accurately flU, 031, 
[[161, and consequently can accommodate arbitrary nonlinear transition surfaces, Lipschitz continuous 
vector fields, and continuous discrete transition maps. We extend this approach using our proposed 
metric topology to prove convergence of simulations to executions of controlled hybrid systems that 



satisfy an orbital stability property described in Section IV Our simulation algorithm may be applied to 
CPS possessing control inputs and overlapping guards, representing a substantial contribution beyond our 
previous efforts [16| and those of others O, 031. 

The remainder of our paper is organized as follows. Section [XT] contains definitions of the topological, 
metric, and dynamical system concepts used throughout the paper. We present our technique for metrization 



and relaxation of controlled hybrid systems in Section III and apply these constructions to define a metric 
for comparing trajectories in controlled hybrid systems. We then develop our algorithm for numerical 
simulation of controlled hybrid system executions in Section |IV} where we apply our trajectory metric 
to prove uniform convergence of simulations to orbitally stable executions. The technique is illustrated 
in Section M using examples for accuracy, verification, and novel controlled hybrid system behavior. 



Section VI contains concluding remarks regarding future research directions. 



II. Preliminaries 
We begin with the definitions and assumptions used throughout the paper. 



A. Topology 

The 2-norm is our finite-dimensional norm of choice unless otherwise specified. Given P, the set of 
all finite partitions of K, and n G N, we define the total variation of / : M. — > M. n by: 

V(/) = su p|E H/(*i+0 - I MT=o G P, m € N J . (1) 

The total variation of / is a semi-norm, i.e. it satisfies the Triangle Inequality, but does not separate 
points. / is of bounded variation if V(/) < oo, and we define BV(El, IR n ) to be the set of all functions 
of bounded variation from K to M 71 . 

Given n e N and D C W 1 , dD is the boundary of D, and int(-D) is the interior of D. Recall that 
given a collection of sets {Sa} ae A> where A might be uncountable, the disjoint union of this collection 
1S LLe.4 $a = [JaeA x a set that is endowed with the piecewise-defined topology. 

The next definition formalizes equivalence relations in topological spaces induced by functions. 

Definition 1. Let S be a topological space, A, B C S two subsets, and f : A — >■ B a function. An 
/-induced equivalence relation is defined as: 

A f ={(a,b) ES xS\ae r\b), or b e r\a), or a = b) . (2) 
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a,b G S are /-related, denoted by a ~ b, if (a, 5) 6 Aj. Moreover, the equivalence class of x G S is 
defined as [x]^ = |a G 5 | a ~ x j. The quotient of S under Af is denoted by and endowed with the 
quotient topology (see Chapter 3 in KTN for more details). 

Note that Af is reflexive, symmetric, and transitive, i.e. an equivalence relation. An important application of 
the function-induced quotient is the construction of a single topological space out of several disconnected 
sets. Indeed, given a collection of sets {S a } a€ A, where A is some index set, and a function f:U—t 

U q& a s <*> where U C H aeA S °" tnen S = s a7 is a topological space. 

Next, we present a useful concept from graph theory that simplifies our ensuing analysis. 

Definition 2. Let (J, Y) be a directed graph, where J is the set of vertices and Y C J x J is the set 
of edges. Then, given j G J, define the neighborhood of j, denoted Afj, by: 

N , j ={eeT\3j'eJs.t.e=(j,f)}. (3) 

B. Length Metrics 

Every metric space has an induced length metric, defined by measuring the length of the shortest curve 
between two points. Throughout this paper, we use induced length metrics to metrize the function-induced 
quotients of disjoint unions of sets. To formalize this approach, we begin by defining the length of a curve 
in a metric space. 

Definition 3. Let (S,d) be a metric space, I C [0,1] be an interval, and 7: / — )■ S be a continuous 
function. Define the length of 7 under the metric d by: 

L d {i) = sup|^d( 7 (t i ),7(t i+1 )) I k G N, {U} k l=0 C /, t < h < ■■ ■ < t k \ ■ (4) 

We now define a generalization of continuous curves for quotiented disjoint unions. 

Definition 4. Let{S a } a< - A be a collection of sets and /:[/—)■ LLe.4 S a , where U C LLe.4 S a . 7 : [0, 1] — >■ 
LLe.4 ^ /-connected if there exists {ti}*L C [0, 1] with = t < t\ < . . . < t k = 1 such that 7|[t ii t l+1 ) 
is continuous for each i G {0, 1, . . . , k — 2}, 7|[t fc _ 1 ,t fc ] is continuous, and lim^ j(t) ~ for each 
i G {0, 1, . . . , k — 1}. Moreover, in that case {U} i=0 is called a partition of 7. 

Note that, since each section j\\t u t i+ i) is continuous, it must necessarily belong to a single set S a for some 
a £ A because the disjoint union is endowed with the piecewise-defined topology. Letting id denote the 
identity function, note that if A contains only one element, then every id-connected curve is simply a 
continuous curve over the original set. Fig. [la] shows an example of a connected curve over a collection 
of two sets. 

Using the concept of connected curves, we now define the induced length distance of a collection of 
metric spaces. The induced length distance is a generalization of the induced metric defined in Chapter 2 

in ma. 



Definition 5. Let {(S a , d a )} aeA be a collection of metric spaces, and let {X a } al - A be a collection of sets 
such that X a C S a for each a G A. Furthermore, let f : U — > LLe.4 X a , where U C LLe^ Xa> an d let 
X = - ■ d { g : X x X — > [0, 00] is the /-induced length distance of X, defined by: 

{k-i 
J2 Ld °i MfcA-n)) I 7: [0, 1] ]J X a , 7(0) = p, 7(1) = q, 
i=0 aeA 

7 is j -connected^ {ti}^ =0 is a partition ofj, (5) 
R}f=o SJ - l\\ti,t i+ i) c X ^ W 
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(a) /-connected curve 7 with partition {U}* =0 , where Si = [0, l] 2 , 82 = [2,3] x (b) Illustration of a controlled hybrid 

[0, 1], and /: {1} x [0, 1] -H>{2} x [0, 1] is defined by ,f(l,x) = f(2,x). system with three modes. 

Fig. 1. An /-connected curve, Fig. [Ta] and an example of a controlled hybrid system, Fig. |lb| 



We use induced length distances in two situations: first, to define a distance over subsets of a single metric 
space, and second, to define a distance over disjoint unions of metric spaces. Definition [5] addresses both 
types of distance functions in a single compact form. Although induced length distances are non-negative, 
symmetric, and subadditive, they do not separate points in general (Section 2.3 in |fT8l ), i.e., every induced 
length distance is a pseudometric, but not necessarily a metric. The following is an important particular 
case of induced length distance. The proof can be found in Section 2.3 in [fT8ll and is omitted here since 
it falls outside the scope of this paper. 

Lemma 6. Let (S, d) be a metric space and X C S. Then d^x, the id-induced length distance on X, is 
a metric. Moreover, the topology on X induced by d^x is equivalent to the topology on X induced by d. 



C. Controlled Hybrid Systems 

Motivated by the definition of hybrid systems presented in [2], we define the class of hybrid systems 
of interest in this paper. 

Definition 7. A controlled hybrid system is a tuple H = {J, T, V, U, J 7 , Q, 71), where: 

• J is a finite set indexing the discrete states of 7i; 

• r C J x J is the set of edges, forming a directed graph structure over J; 

• V ={Dj}. e j is the set of domains, where each Dj is a subset ofWL nj , rij 6 N; 

• U C lR m is the range space of control inputs, m e N; 

• T = {fj}j e j i s tn e set of vector fields, where each fj-.WLxDjXU — > M." j is a vector field defined 
on Dj; 

• G = {Ge} e£ r zs tne set of guards, where each Gq,/) C dDj is a guard in mode j E J that defines 
a transition to mode j' G J ; and, 

• 1Z = {i? e } egr is the set of reset maps, where each R(jj') : G(j,f) ~ > Dj' defines the transition from 
guard Gijjiy 

For convenience, we sometimes refer to controlled hybrid systems as just hybrid systems, and we refer 
to the distinct vertices within the graph structure associated with a controlled hybrid system as modes. 
Note that each domain in the definition of a controlled hybrid system is a metric space with the Euclidean 



distance metric. A three-mode controlled hybrid system is illustrated in Fig. lb 
Next, we impose several technical assumptions. 

Assumption 8. For each j G J, fj is Lipschitz continuous. That is, there exists L > such that for each 
j G J, ti,t 2 G R, Xi,x 2 G Dj, and Ui,u 2 G U: 

||/ i (ti,xi,'Ui) - fj(t2,X2,u 2 )\\ < - t 2 \ + ||aei - x 2 \\ + |K - u 2 \\). (6) 
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Assumption [8] guarantees the existence and uniqueness for ordinary differential equations in individual 
domains. A proof of this result can be found in Section 2.4.1 in [fT9ll . 

Lemma 9. Let H be a controlled hybrid system. Then for each j G J , each initial condition p G Dj, and 
each control u G BV(M., U), there exists an interval / CM with G J such that the following differential 
equation has a unique solution: 

x(t) = fj(t,x(t),u(t)), tel, x{0)=p. (7) 

x is the integral curve of fj with initial condition p and control u. 

The following definition is used to construct executions of a controlled hybrid system. 

Definition 10. Let H be a controlled hybrid system, j G J , p G Dj, and u G BV(J5&, U). x: I — » Dj is 
the maximal integral curve of fj with initial condition p and control u if given any other integral curve 
with initial condition p and control u, such as x: I — > Dj, then I C I. 

Assumption 11. Let % be a controlled hybrid system. Then the following statements are true: 

(1) For each j G J , Dj is a compact, rij-dimensional, manifold with boundary. 

(2) U is compact. 

(3) For each e G V, G e is a closed, embedded, codimension 1, submanifold with boundary. 

(4) For each e G T, R e is continuous. 



Given a maximal integral curve x: I — > Dj, a direct consequence of Definition 10 and Assumption 11 



is that either sup / = +oo, or sup / = t' and x(t') G dDj. This fact is critical during the definition of 



executions of a controlled hybrid systems in Section IV 



III. Metrization and Relaxation of Controlled Hybrid Systems 

In this section, we metrize a unified family of spaces containing all the domains of a controlled hybrid 
system H. The constructed metric space has three appealing properties: 

(1) the distance between a point in a guard and its image via its respective reset map is zero; 

(2) the distance between points in different domains are properly defined and finite; and, 

(3) the distance between points is based on the Euclidean distance metric from each domain. 



A. Hybrid Quotient Space 

Using Definitions [5] and [7J we construct a metric space where the executions of a controlled hybrid 
system reside. The result is a metrization of the hybrifold fl2). 

Definition 12. Let 1-ibe a controlled hybrid system, and let 

R: ]JGW jjDj (8) 

eer jej 

be defined by R(p) = R e (p) far each p G G e . Then the hybrid quotient space ofH is: 

M = Hq^. (9) 



Fig. [2] illustrates the details about the construction described in Definition 12 
The induced length distance on M. is in fact a distance metric: 

Theorem 13. Let H be a controlled hybrid system, and let di t M be the R-induced length distance of A4, 
where R is defined in ([8]). Then d^M is a metric on Ai, and the topology it induces is equivalent to the 
R-induced quotient topology. 
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Fig. 2. The disjoint union of D\ and D2 (left) and the hybrid quotient space M obtained from the relation (right). 



Proof: We provide the main arguments of the proof, omitting the details in the interest of brevity. First, 
note that each domain is a normal space, i.e. every pair of disjoint closed sets have disjoint neighborhoods. 
Second, note that each reset map is a closed map, i.e. the image of closed sets under the reset map are 
closed. This fact follows by Condition [3] in Assumption 1 1 since each guard is compact, thus reset maps 



are closed by the Closed Map Lemma (Lemma A. 19 in II2010 . 

Let D = Yijej ®j an d p,q E D. Then their equivalence classes, [p]^ and [q]^, are each a finite collection 
of closed sets. Moreover, since we can construct disjoint neighborhoods around each of these closed sets, 
then we can conclude that there exists 5 > such that di,M ([p]^ ,[<?]#) > The proof concludes by 
following the argument in Exercise 3.1.14 in lfT8l . i.e. since each connected component in D is bounded, 
then Ai is also bounded (in the quotient topology). Then, by a straightforward extension of Theorem 5.8 
in IfTTTl . we get that the identity map from Ai to the space constructed by taking the quotient of all the 
points in D such that i» has zero distance is a homeomorphism, thus they have the same topology. 
Note that, by using the metric g^x, all i?-connected curves are continuous. As we show in Section 



IV 



this implies that executions of controlled hybrid systems are continuous in the topology induced by d^M- 

B. Relaxation of a Controlled Hybrid System 

To construct a numerical simulation scheme which does not require the exact computation of the time 
instant when an execution intersects a guard, we require a method capable of introducing some slackness 
within the computation. This is accomplished by relaxing each domain along its guard and then relaxing 
each vector field and reset map accordingly in order to define a relaxation of a controlled hybrid system. 

To formalize this approach, we begin by defining the relaxation of each domain of a controlled hybrid 
system, which is accomplished by first attaching an £-sized strip to each guard. 

Definition 14. Let % be a controlled hybrid system. For each e E V, let S e e = G e x [0,e] be the strip 
associated to guard G e . For each j E J, let 



]\ G e -+ ]\ SI (10) 

e&A/",- e&ATi 



be the canonical identification of each point in a guard with its corresponding strip defined for each 
p G G e as Lj{p) = {jp, 0) G Si. Then, the relaxation of Dj is defined by: 



^Ii(lI e6Afj Sf) 



D- = (11) 



By Condition in Assumption 1 1 each point on a strip SI of Dj is defined using rij coordinates 



(Ci, • • • , Cn — i) T), shortened (£, rJTwhere r is called the transverse coordinate and is the distance along 



the interval [0, e]. An illustration of Definition 14 together with the coordinates on each strip is shown in 
Fig. [3a 



(a) Disjoint union of D\ and the strips in its neighborhood, {5 , |} egA r i (left), and the 
relaxed domain D\ obtained from the relation A tl (right). 



(b) Relaxed vector field ff on Df . 



Fig. 3. An illustration of the construction of a relaxed domain, Fig. [3a| and the relaxed vector field defined on it, Fig. |3b| 



We endow each SI with a distance metric in order to define an induced length metric on a relaxed 
domain D £ . 

Definition 15. Let j G J and e G A/}. Endow Dj with d^Dj as its metric, and : S e £ x SI —> [0, oo) as 
the metric on S e e , defined by: 

ds i ((C,T),(C,r'))=d hGe (CX') + \r-r'\. (12) 

We now define a length metric on relaxed domains using Definitions [4] and [5J 

Theorem 16. Let j G J , and let di £><? be the L~-induced length distance on D% where ij is as defined 
in ( |10] ). Then di^ is a metric on D e , and the topology it induces is equivalent to the ij-induced quotient 
topology. 

Proof: Since d ijD e is non-negative, symmetric, and subadditive, it remains to show that it separates 
points. Let p,q G D £ . First, we want to show that [p] t . = [q] L , whenever d iiD e(p,q) = 0. Note that for 
each e E J\fj and each pair p,q G G e , and by the Definition [f>j and 15 ds^(p,q) > di t D 3 (p,q), hence no 



connected curve that "jumps" to a strip can be shorter than a curve that stays in Dj. This fact immediately 
shows that for p,q<E Dj, d ijD ?(p, <?) = implies that [p] t . = [q] L .. The case when one of the points is in 
G e x (0,e] C Si follows easily by noting that those points can be separated by a suitably-sized ds^-ball. 
The proof concludes by following the argument in Exercise 3.1.14 in [fT8ll . as we did in the proof of 
Theorem [T3l ■ 
Refer to di D e as the relaxed domain metric. Note that Theorem 



as in the proof of Theorem 13 but we prove Theorem 16 to emphasize the utility of the inequality relating 



16 can be proved using the same argument 



the induced metric on a domain and the metric on each strip. 
Next, we define a vector field over each relaxed domain. 

Definition 17. Let j G J. For each e G A/}, let the vector field on the strip S e e , denoted f e , be the unit vector 
pointing outward along the transverse coordinate. In coordinates, f e (t, ((,t),u) = (0, 0, l) .Then, 

f coords. 



the relaxation of /, is: 



fj(t,x,u) ifxeDj, 



j/ ? f e (t,x,u) if x G G e x (0,e] C S £ , for some e G A/}. ^ ^ 

Note that the relaxation of the vector field is generally not continuous along each G e , for e G A/}. As 
we show in the algorithm in Fig. [7J this discontinuous vector field does not lead to sliding modes on the 
guards 11211 . [|22l . since the vector field on the strips always points away from the guard. An illustration 



of the relaxed vector field /J on D e is shown in Fig. 3b 
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Fig. 4. The disjoint union of D\ and _D| (left), and the relaxed hybrid quotient space A4 e obtained from the relation A^ e (right). 



The definitions of relaxed domains and relaxed vector fields allow us to construct a relaxation of the 
controlled hybrid systems as follows: 

Definition 18. Let % be a controlled hybrid system. We say that the relaxation of H is a tuple 'W = 
{J, T, V £ , U, T e , G £ , K E ), where: 

• V s = is the set of relaxations of the domains in V, and each D| is endowed with its 
induced length distance metric d^op' 

• T e = \fj}j e j is tne set of relaxations of the vector fields in T; 

• Q £ = {G e e } e€T is the set of relaxations of the guards in Q, where G £ e = G e x {e} C S £ for each 
e G T; and, 

• lZ e = {-Rg} e r i s the set of relaxations of the reset maps in 1Z, where R e e : G e e — > Djr for each 
e=(j,f)er and Rl((,e) = R e (() for each (eGl. 



C. Relaxed Hybrid Quotient Space 



Analogous to the construction of the metric quotient space Ai, using Definitions [5] and 18 



we construct 

a unified metric space where executions of relaxations of controlled hybrid systems reside. The result is a 
metrization of the hybrid colimit [3] rather than a metrization of the hybrifold as in the previous section. 

Definition 19. Let H £ be the relaxation of the controlled hybrid system H. Also, let 

R e \\Gl^\\D] (14) 

eer jej 

be defined by R £ (p) = R e e (p) for each p G G £ e . Then the relaxed hybrid quotient space of W is: 

AT = %^5. (15) 

% 

The construction in Definition 19 is illustrated in Fig. |4j 

We now show that the induced length distance on A4 e is indeed a metric. We omit this proof since it 
is identical to the proof of Theorem [T3j 

Theorem 20. Let 7-L be a controlled hybrid system, let 7-L £ be its relaxation, and let d it M e & e the R e - 
induced length distance of Ai e , where R £ is defined in (14). Then di t M E is a metric on A4 e , and the 
topology it induces is equivalent to the Rr-induced quotient topology. 

All _R e -connected curves are continuous under the metric topology induced by d{ ^ which is important 
when we study executions of hybrid systems in Section |IV} 

As expected, the metric on M. e converges pointwise to the metric on Ai. 
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Theorem 21. Let % be a controlled hybrid system, and let V. £ be its relaxation. Then for all p,q G M., 

di,M s (Pi <?) -> d i)M (jP, q) as e -)■ 0. 

Proof: Abusing notation, let 1/(7) denote the length of any connected curve 7, defined as the sum of 
the lengths of each of its continuous sections under the appropriate metric. First, note that d^MiPiO) < 



di t M e {PiQ)- This inequality follows since, as we argued in the proof of Theorem 16 given an edge 



(j)j') e T, dge , t (p',q') > d iiDj (p',q') for any pair of points p',q' G G(jj>y Thus, adding the strips 
{S e } ee in Ai £ only make the length of a connected curve longer. 

Now let D = Yijej ®3 an( ^ ® £ = Uj & jDj. Given 5 > 0, there exists 7: [0, 1] — > D, an i?-connected 
curve with partition {U} k i=0 , such that 7(0) = p, 7(1) = q, and d^MiPiQ) < -^(7) < d itM (p,q) + 5. 
Moreover, without loss of generality let j £ : [0, 1] — >■ _D e be an i? e -connected curve that agrees with 7 on 
_D, i.e. each section of 7 £ on 2) is identical, up to time scaling, to a section of 7 Thus 7 s has at most k 
^-length extra sections, and £(7) < L(Y) < £(7) + Thus, d;,.M e (p, 9) < ^(7 £ ) < ^t,x (p> 9) + + 5- 
The result follows after noting this inequality is valid for all 8 > 0, thus d^M^ip, q) < di^ip, <?)• ■ 



Note that Theorem 21 does not imply that the topology of M. £ converges to the topology of M.. On the 
contrary, A4 £ is homotopically equivalent to the graph ( J, V) for each e > 0, whereas the topology 
of M. may be different 0. 

We conclude this section by introducing metrics between curves on Ai £ . 

Definition 22. Let I C [0, 00) a bounded interval. Given any two curves 7,7': I — > Ai £ , we define: 

pf(7,y) = snp{d iM e (7(^,7' (t)) I t g /}. (16) 



Our choice of the supremum among point-wise distances in Definition [22] is inspired by the sup-norm 
for continuous real-valued functions. 

IV. Relaxed Executions and Discrete Approximations 

This section contains our main result: discrete approximations of trajectories of controlled hybrid 
systems, constructed using any variable step size numerical integration algorithm, converge uniformly 
to the actual trajectories. This section is divided into three parts. First, we define a pair of algorithms that 
construct executions of controlled hybrid systems and their relaxations, respectively. Next, we develop 
a discrete approximation scheme for executions of relaxations of controlled hybrid systems. Finally, we 
prove that these discrete approximations converge to the executions of the original, non-relaxed, controlled 



hybrid system using the metric topologies developed in Section III 



A. Execution of a Hybrid System 

We begin by defining an execution of a controlled hybrid system. This definition agrees with the 
traditional intuition about executions of controlled hybrid systems which describes an execution as evolving 
as a standard control system until a guard is reached, at which point a discrete transition occurs to a new 
domain using a reset map. We provide an explicit definition to clarify technical details required in the 
proofs below. Given a controlled hybrid system, H, as in Definition [7J the algorithm in Fig. [5] defines 
an execution of % via construction. A resulting execution, denoted x, is an /^-connected curve from 
some interval / C [0, 00) to Yijej Thus, abusing notation, we regard x as a continuous curve on M. 



Abusing notation again, we regard 1 as a piece-wise continuous curve on A4 £ for each e > 0. Fig. 6a 
shows an execution undergoing a discrete transition. 

Note that executions constructed using the algorithm in Fig. [5] are not necessarily unique. Indeed, 



Definition 10 implies that once a discrete transition has been performed, the execution is unique until a 
new transition is made; however, the choice in Step [8] is not necessarily unique if the maximal integral 
curve passes through the intersection of multiple guards. It is not hard to prove that a sufficient condition 
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Require: t = 0, j G J, p G D jt and u G BV(R, U). 
i: Set sc(0) = p. 

2: loop 

3: Let 7: 7 — > £)j be the maximal integral curve of fj with control u such that j(t) = x(t). 

4: Let t' = sup / and x(s) = 7(5) for each s G [t,t'). > Afote j/i' < 00, 7(f) G ftDj-. 

5: if f = 00, or G A/} such that 7(f) G G e then 

6: Stop. 

7: end if 

8: Let (j,f) G A/j be such that j(t') G Gyj'). 

9: Set x{t') = R (j:f) (7(f)), * = and J = > Note j(t') ~ 

10: end loop 

Fig. 5. Algorithm to construct an execution of a controlled hybrid system H. 




for uniqueness of executions is that all the guards are disjoint, even though, as we show in Section V-C 
uniqueness of the executions can be obtained for some cases where guards do intersect. 

With the definition of execution of a controlled hybrid system, we can define a class of executions 
unique to controlled hybrid systems. 

Definition 23. An execution is Zeno when it undergoes an infinite number of discrete transitions in a 
finite amount of time. Hence, there exists T > 0, called the Zeno Time, such that the execution is only 
defined on I = [0, T). 

Zeno executions are hard to simulate since they apparently require an infinite number of reset map 
evaluations, an impossible task to implement on a digital computer. A consequence of the algorithm in 
Fig. [5] is that if x: I — > M. is an execution such that T = sup / < 00, then either 

(1) x has a finite number of discrete transitions on / = [0, T], and x(T) G dDj for some j G J, or 

(2) a; is a Zeno execution and / = [0, T). 

We now introduce a property of Zeno executions of particular interest in this paper: 

Definition 24. Let H be a controlled hybrid system, p G M, u G BV(M, U), and x: [0, T) — > M be 
a Zeno execution with initial condition p, control u, and Zeno Time T. x accumulates at p' G M. if 
lim t ^ T d iM (x(t),p') =0. 



Examples of Zeno executions that do not accumulate can be found in [23|. Fig. 6b shows a Zeno execution 
that accumulates at p'. Note that for p' to be a Zeno accumulation point, it must belong to a guard of a 
controlled hybrid system. 

Since Ai is a metric space, we can introduce the concept of continuity of a hybrid execution with 
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Require: t = 0, j G J, p G D jt and u G BV(R, U). 

i: Set x e (0) = p. 

2: loop 

3: Let 7: 7 the maximal integral curve of with control u such that j(t) = x £ (t). 

4: Let t' = sup / and x £ (s) = 7(3) for each s G [£,£')• ■> < °°> men e 

5: if f = 00, or $e G Mj such that j(t') G G e then 

6: Stop. 

7: end if 

8: Let (j,f) G Mj such that 7(f) G G ijif) , thus (7(f), 0) G 

9: Set x e (t'+r) = (7(f), r) for each r G [0, e). 

10: Set x £ (t'+e) = . /) (7(t'),e), t = t'+e, and j = j'. > Note (7(f), £ ) ~V(t'+e). 

ii: end loop 

Fig. 7. Algorithm to construct a relaxed execution of a relaxation of a controlled hybrid system, H £ . 



respect to its initial condition and control input in a straightforward way. Employing this definition, we 
can define the class of executions that are numerically approximable: 

Definition 25. Let H be a controlled hybrid system, and assume that all the executions of H are unique. 
Denote by X( p n ) : I( PtU ) «M the hybrid execution of H with initial condition p G M. and control 
u G BV{R,U). Given T > 0, we say that the map (p,u) 1— > X(p, u ) is orbitally stable in [0, T] at 
{p',u r ) G M x BV(R,U) if there exists a neighborhood of (p',u'), say AVyx) x BV(R,U), such 

that the following conditions are satisfied: 

(1) [0,T] C I( p , u )for each (p,u) G N {p , y) . 

(2) The map (p,u) i-> xu, )U )(t) is continuous at (p',v!) for each t G [0, T]. 

As observed by [|24l|. orbitally stable executions are exactly the type of executions of a controlled hybrid 
system that can be approximated. Indeed, if an execution is not orbitally stable then there exists a time 
f such that, if another execution is initialized arbitrarily close to x(t') or if an arbitrarily close control 
is applied, then the pair of executions travel through different sequences of discrete transitions and can 
never be made arbitrarily close. Fig. 6c shows a non-orbitally stable execution that intersects the guard 
tangentially. 



B. Relaxed Execution of a Hybrid System 

Next, we define the concept of relaxed execution for a relaxation of a controlled hybrid system. The 
main idea is that, once a relaxed execution reaches a guard, we continue integrating over the strip with 
the relaxed vector field, f e , as in Definition 17 Given the controlled hybrid system, H, its relaxation, H e , 
for some e > 0, the algorithm in Fig. [7] defines a relaxed execution of H 6 via construction. The resulting 
relaxed execution, denoted x £ , is a continuous function defined from an interval / C [0, 00) to Ai £ . Note 
that this algorithm is only defined for initial conditions belonging to Dj for some j G J since the strips 
are artificial objects that do not appear in H. The generalization to all initial conditions is straightforward; 
we omit it to simplify the presentation. 

Step [9] of the algorithm in Fig. [7] relaxes each instantaneous discrete transition by integrating over the 
vector field on a strip, hence forming a continuous curve on Ai £ . Also note that our definition for the 
relaxed execution over each strip Si, also in Step |9| is exactly equal to the maximal integral curve of f e . 
Fig. 8a shows an example of a relaxed mode transition produced the algorithm in Fig. [7J Given a hybrid 
system H and its relaxation T-L £ , the relaxed execution of 1-L £ produced by the algorithm in Fig. [7] is a 
delayed version of the execution of % produced by the algorithm in Fig. |5J since the relaxed version has 
to expend e time units during each discrete transition. In that sense, our definition of relaxed execution 
is equivalent to an execution of a regularized hybrid systems 0. 
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Note that if a relaxed execution is unique for a given initial condition and input, then the corresponding 
hybrid execution is also unique, but not vice versa. Indeed, consider the case of a hybrid execution 
performing a single discrete transition at a point, say p, where two guards intersect, i.e. p E G e and 
p E G e i, such that R e (p) = R e '{p)- In this case the hybrid execution is unique, but its relaxed counterpart 
either evolves via S e or S e i, hence obtaining 2 different executions. Nevertheless, both relaxed executions 
reach the same point after evolving over the strip. 

Next, we state our first convergence theorem. 

Theorem 26. Let 7-Lbe a controlled hybrid system, and H, £ be its relaxation. Let p E Ai, u E BV(JRL, U), 
x : I — > Ai £ be an execution of H with initial condition p and control u, and let x £ : I s — > Ai £ be a 
corresponding relaxed execution of x. Assume that the following conditions are satisfied: 

(1) x is orbitally stable with initial condition p and control u; 

(2) x has a finite number of discrete transitions or is a Zeno execution that accumulates; and 

(3) there exists T > such that for each e small enough, [0, T] C I C\ I £ if x has a finite number of 
discrete transitions, and [0, T) C / fl I £ if x is Zeno. 

Then, lim e ^ pf 0)T ] {x,x £ ) =0. 

Proof: We provide the main arguments of the proof, omitting some details in the interest of brevity. 
First, given j E J and [r, r') C [0, T] such that x(t) E Dj for each t E [r, r'), then, since x|[ r>r /) is 
absolutely continuous, for each t,t' E [t, t'), 

d iM e[x{t),x{t')) < L diDj (x\[ tit ')) = J \\fj(s,x(s),u(s))\\ ds < K(t' - t), (17) 

where K = sup {\\f £ (t,x,u)\\ \ j Ej, tE [0,T], x E M £ ,u E U) < oo. 

Second, let k E N and {Aj} 4 fe =0 C [0, 1] be a sequence such that = A < Ai < . . . < \ k = 1. Given 



e > 0, let 7 t : [0,1] — > Ai £ be defined by 7t(A) = x Xe (t). Thus, by Theorem 21 and the algorithm in 
Figure 7} 7 t (0) = x°(t) = x(t) and 7t(l) = x £ {t). Assume that x £ (t) E Dj for each t E [r + e,r' + e), 
where [r, r') is as defined above. Using Picard's Lemma (Lemma 5.6.3 in ll25ll ). for each t E [r + e, t'), 

\\x £ (t + e) - x(t) || < e L[t - T) (\\x £ (t + e) - x(r) \\ + 

+ J \\fj(s,x(s),u(s)) — fj(s + e,x(s),u(s + ds 

< e L{t - T) (\\x £ (t + e) - x(r) \\+L j^e+ \\u{s) - u{s + e) \\ ds^j 

< e L ^ (\\x £ (t + e) - x{t) \\ + (L + V(w)) {t - r)e) , 



(18) 
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where we have used a standard property of the functions of bounded variation (Exercise 5.1 in Il26l0 . Thus, 
if we assume that \\x e {r + e) — x(t)\\ = 0(e), i.e. that there exists C > such that ||a; e (r + e) — x(t)\\ < 
Ce, then \\x e (t + e) — x(t)\\ = 0(e) for each t G [r + e,r'). Using the same argument as above 
i^Ai+ie^ + e) — a; Ai£ (t)|| = 0((\ i+ i — \i)e), which implies that is continuous for each t G [r + e, t'), 
and that L(7 t ) = 0(e), hence d itD] (x e (t + e),x(t)) = 0(e). 

Assuming now that x performs 2 discrete transitions at times r, r' G [0,T], such that r + e < r', 
transitioning from mode j to j', and the from mode f to j". Note that, by definition, x|r 0)T ) = x E |[ 0>r ). 
Moreover, since £ is orbitally stable, we know that x £ performs the same 2 discrete transitions for e small 
enough. Let r £ + e G [0,T] be such that x e (r e + e) G G(j',j")- Note that |r E — r'| = 0(e) since x e — > x 
uniformly and x is Lipschitz continuous (both propositions shown above). Assume that r' < r £ + e and 
consider the following upper bounds: 

(1) lfte[r,r + e), then x(t) G D, y and x e (t) G thus: 

<W (z(t), x £ (t)) < x(r)) + ds« _ , /} (x(r), x £ (t)) = 0(e) . (19) 

(2) If t G [r + e, r'), then x(t), x £ (t) G Dy, thus, using the bound obtained above: 

d iM e (x(t),x £ (t)) < d i>D ., (x(t),x(t - e)) + d ijDjl (x(t - e),x £ (t)) = 0(e) . (20) 

(3) If t G [r', t £ + e), then x(t) G Dyi and x £ G Dy, thus, denoting lim^' x(t) = x(r'~): 

d iM e(x(t),x £ (t)) < d i>D (x(t),x(r')) + d« . (x(r'), x(r /_ )) + 

+ (x(t'-),x £ (t £ + e)) + d i>Djl (x £ (t £ + e),x £ (t)) (21) 
< 0(e) . 

(4) If t G [r £ + e, r £ + 2e), then x(t) G -Dj» and x e G S'L.y/), thus: 

di^- (a:(t), x £ (t)) < d i>D .„ (x(t),x(r')) + dsy, ^ {x(t') , x £ (t)) < 0(e) . (22) 

(5) If t G [t £ + 2e,T], then x(t) 1 x £ (t) G -Dj", thus we get the same bound as in case[5J 
Therefore, p £ T j(x, x £ ) = 0(e) as desired. Note that the general case, with an arbitrary number of discrete 
transitions and where the discrete transitions of x £ occur before the discrete transitions of x, follows by 
using the a similar argument as above by properly considering the time intervals and then applying the 
upper bounds inductively. 

Next, let us consider the case when x is a Zeno execution that accumulates on p'. Let 5 > 0, then 
%\{o,t-5] has a finite number of discrete transitions, and as shown above, d it M e (x(T—S),x £ (T—5)) = 0(e). 
Moreover, d iM s(x(T- 6), x(t)) =0(5) and d iM ,(x £ (T - 5),x £ (t)) =0(5) for each t G [T-5,T).The 
conclusion follows by noting that these bounds are valid for each 5 > 0. ■ 

C. Discrete Approximations 

Finally, we are able to define the discrete approximation of a relaxed execution, which is constructed 
as an extension of any existing ODE numerical integration algorithm. Given a controlled hybrid system 
"H, Aj'. E x ]R n J x U — > IR nj , where h > and j 6 J, is a numerical integrator of order to, if given 
p G Dj, u G BV(R,U), x the maximal integral curve of fj with initial condition p and control u, 
N = |_^J, and a sequence {z k }^ =0 with z = p and Zk+\ = A^{kh, Zk,u(kh)), then sup{ \\x(kh) — z k \\ 
k G {0, . . . ,iV}} = 0(lf). This definition of numerical integrator is compatible with commonly used 
algorithms, including Forward and Backward Euler algorithms and the family of Runge-Kutta algorithms 
(Chapter 7 in 112710 . The algorithm in Fig. [9] defines a discrete approximation of a relaxed execution of 
H 6 . The resulting discrete approximation, for a step size h > 0, denoted by z e,h , is a function from a 
closed interval / C [0, oo) to M e . 
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Require: h > 0, k = 0, j G J, and p <E Dj. 

l: Set t = and 2 E ' h (0) = p. 

2: loop 

3: Find n' = mf{n G N | Af~" (t k , z £ > h {t k ),u{t k )) G D £ ) 

4: if n' = oo then 
5: return z £ > h \ [0ttk ]. 

6: end if 

7: Set £ fe+1 = t fc + h2~ n ', and ^' h (t fc+ i) = *4f ""' (t k , z £ ' h (t k ),u(t k )) . 

8: Set z £ ' h {t) = £+£- z ^(t k ) + -*=^^."(t fc+1 ) for each t G [t fe , t k+1 }. 

L k-\-l L k L k-\-l L k 

9: if 3(j,f) G A/} such that z £ > h {t k+1 ) G S^.,) then 

10: Set (q, r) = z s ' fc (tjb + i) G S^,), and t fc+2 = t k+1 + e - r. 

ll: Set z £>/l (t) = (q,t- t k+1 + r) for each t G [^+1,^+2). 

12: Set z £ - /l (t fc+2 )=i?^, ) (g, £ ), k = k + 2, and j = f. > Note (q, e) ~ z £ ' h (t k+2 ). 

13: else 

14: Set k = k + 1. 

15: end if 

16: end loop 

Fig. 9. Discrete approximation of a relaxed execution of the relaxation of a controlled hybrid system H e . 



We now make several remarks about the algorithm in Fig. |9} First, the condition in Step [4] can only be 
satisfied, i.e. the Algorithm only stops, if z £ ' h (t k ) G dDj and fj(t k ,z £ ' h (t k ),u(t k )) is outward-pointing, 
since otherwise a smaller step-size would produce a valid point. Second, the function z £,h is continuous 
on M. £ . Third, and most importantly, similar to the algorithm in Fig. [7j the curve assigned to z £,h in 
Step 11 is exactly the maximal integral curve of f e while on the strip. By relaxing the guards using strips, 
and then endowing the strips with a trivial vector field, we avoid having to find the exact point where the 
trajectory intersects a guard. Our relaxation does introduce an error in the approximation, but as we show 



in Theorem |27} the error is of order e. Fig. [8b] shows a discrete approximation produced by the algorithm 
in Fig. [9] as it performs a mode transition. 

Theorem 27. Let %be a controlled hybrid system and 1-L £ its relaxation. Let p G M., u G BV(WL, U), and 
let x: I Ai £ be a orbitally stable execution ofH with initial condition p and control u. Furthermore, 
let x £ : I s — > Ai £ be a relaxed execution with initial condition p and control u, and let z £,h : I £,h — > Ai £ 
be its discrete approximation. If [0, T] C I s R I £ ' h for each e and h small enough, then there exists C > 
such that lim^o pf j>] (x £ , z £ ' h ) < Ce. 

Proof: As we have done with the previous proofs, we only provide a sketch of the argument in the 
interest of brevity. Assume that x £ performs a single discrete transition in the interval [0,T] for each e 
small enough, crossing the guard Gu^\ at time r £ . Then, since x is orbitally stable and A h is convergent 
with order cu, for e and h small enough z £ ' h also crosses guard Guji\ at time t^J 1 G [t k /,t k i + i) for some 
k' G N, where {t k }^ =0 is the set of time samples associated to z £)h . Moreover, since x £ (fS) = z £ ' h (0), then 
for each 5 > 0, |r E - t k , +1 \ <5 + 0(h") and \t k , +2 - t £ + e\ = 0(h"), 
Define the following times: 

a m = mm{t k , +1 ,T £ }, a M = max{t k > +1 , r £ }, 
v m = mm{t k , +2 ,T £ + e}, v m = max{t fc / +2 , r e + e}, 

and, in order to simplify our argument, assume that <tm < v m . Then on the interval [0, er m ) we get 
convergence due to A h . On the interval [a m , <ta/) one execution has transitioned into a strip, while the 
other is still governed by the vector field on Dj. On the interval [a M ,cu m ) both executions are inside the 
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strip, and on the interval [cu m , lom) one execution has transitioned to a new domain, while the second is 
still on the strip. After time com both executions are in a new domain, and we can repeat the process. 
Consider the following cases: 

(1) By the convergence of algorithm A h , d i:M ^x £ (a m ), z £ ' h (a m )) = 0{h"). 

(2) Using ([17]) from the proof of Theorem [26} 



d iM ?(x £ ((TM),Z £ ' h ((TM)) < d iM e(x e (a M ),X e {(?m)) + ((T m ) , Z £ ' k ((T m )) + 

+ d iiM e(z E ' h (a m ), z £ ' h (a M )) 

(3) Using the same argument as in the inequality above, 

d iM e (x £ {u m ), z £ ' h (u m )) < d iM e (x £ {a M ), z £ ' h (a M )) + 2e. 

(4) Finally, again using the same argument as in case [2} 

d iM e(x £ {is M ),z e ' h {v M )) < d iM e(x £ (v m ),z £ > h (v m )) +0(h"). 



0{h u ). (24) 



(25) 



(26) 



The generalization to any relaxed execution defined on A4 £ and its discrete approximation follows by 
noting that they perform a finite number of discrete jumps on any bounded interval and that 5' can be 
chosen arbitrarily small. ■ 



Next, we state the main result of this Section, which is a result of Theorems 26 and 27 



Corollary 28. Let %be a hybrid dynamical system and T-L £ be its relaxation. Let p G M., u G BV{R, U), 
x: I — > Ai £ be an execution ofH with initial condition p and control u, x £ : I £ — >■ Ai £ be its corresponding 
relaxed execution, and z £,h : I £,h — > A4 £ be its corresponding discrete approximation. If the following 
conditions are satisfied: 

(1) x has a finite number of mode transitions or is a Zeno execution that accumulates, 

(2) x is orbitally stable, 

(3) [0, T] C / fl I £ H I £ ' h for each e and h small enough, 
then lim e _>o pf n n (x, z £,h ) = 0. 



Proof: Note that, by Theorem 26 together with the Triangle Inequality, this corollary is equivalent 



to proving that p £ j(x £ , z £ ' h ) — > as both e, h — >• 0. Hence we show that p £ j(x £ , z £ ' h ) converges uniformly 
on h as e — >■ 0. Using an argument similar to the one in the proof of Theorem 7.9 in [28], proving the 
uniform convergence on h is equivalent to showing that lim^o limsup^Q p\ (x e , z £ ' h ) = 0, but this is 
clearly true by Theorem [271 therefore obtaining our desired result. ■ 



V. Examples 

We present three examples where our discrete approximation algorithm is applied, first comparing its 



performance to existing state-of-the-art algorithms in Sec. V-A, then showing its use while computing 



the trajectories of a legged-locomotion model in Sec. V-C, and finally implementing a benchmark hybrid 



system verification example in Sec. V-B 



A. Forced Linear Oscillator with Stop 

We consider a single degree-of-freedom oscillator consisting of a mass that is externally forced and 



can impact a plane fixed rigid stop, as in Fig. 10a The state of the oscillator is the position, x(t) G K, 
and velocity, x(t) G K, of the mass. The oscillator is forced with a control u G BV(M, R). The oscillator 
is modeled as a controlled hybrid system with a single mode, denoted D, and single guard corresponding 
to the mass impacting the stop with non-negative velocity, denoted G: 

D = { (x(t),x(t)) G R 2 | x(t) < x max } 
G = { (x(t),x(t)) G M 2 | x(t) = x max , x(t) > 0} 



(27) 
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m 



x(t) 




(a) Forced linear oscillator with stop. (b) Position of the analytical solution of Ex- (c) Position of the analytical solution of Ex- 

ample 1 in Table [I] (solid line), and posi- ample 2 in Table [I] (solid line), and posi- 
tion of the stop (dotted line). tion of the stop (dotted line). 




(d) Comparison of the algorithm in Fig. |9] vs. (e) Comparison of the algorithm in Fig. [9j vs. (f) Computation times of the algorithm in 
the PS Method for the examples in Table|I] the PS Method for the examples in Table]!] Fig. [9] vs. the PS Method for the examples 
using pfo,^] ■ using p. in Table [I] 

Fig. 10. A mechanical system (Fig. |10a} an d a pair of examples (Figs. [T0b"] an d|10cfr chosen to illustrate the accuracy of the algorithm in 
Fig. [9] vs. the PS Method (Figs. |10d| and |10e[ l and their computation times (Fig. |10f[ >. 



Upon impact, the state is updated using the reset map R(x,x) = (x,—cx), where c G [0,1] is the 
coefficient of restitution. Within the single domain, the dynamics of the system are governed by x(t) + 
2ax(t) + oo 2 x(t) = m -1 w(t), where u = y/mr^k, a = 0.5 m' 1 jj,, k is the spring constant, and /i is the 
damping coefficient. 

Given an initial condition (x(t ), x(t )) = (xo,Xo) G D, the oscillator's motion is analytically de- 
termined by x(t) = e~ at (A n cos(Cot) + B n sm(Cot)) + Co" 1 L u(s)e~ a ^~ s ^ sm(p(t — s)) ds for each t G 
[t n _i, t n ), where Co = y/to 2 — a 2 (assuming that the damping is sub-critical), with t n such that x(t~) = x max 
for each n G N, and A n and B n are determined by the given initial conditions when n = 0, or those 
determined by applying the reset map to x(t~) when n > 1. Note that determining the impact times can 
be done analytically. The analytical solution holds provided that the mass does not stick to the stop, since 
in that case the dynamics are given by x(t) + 2ax(t) + to 2 x(t) = m" 1 (u(t) + \(t)), where A(t) G K 
denotes the force generated by the stop to prevent movement. This equation holds as long as x(t) = x max , 
x(t) = x(t) = 0, and the reaction of the stop is negative, i.e. A(t) > mco 2 x max . For the contact to cease, 
A(t) — mco 2 x max must become zero and change sign. Once this happens, the analytical solution can be 
used again to construct the motion of the mass with the initial condition (x max ,0). 

Assuming that the forcing u is continuous (an assumption that is violated by many control schemes 
such as ones generated via optimal control) a convergent numerical simulation scheme, which we call the 
PS Method, to determine the position of a mechanical system with unilateral constraints was proposed 
in lfl4ll . Fixing a step-size h > 0, their approach is a two-step method that for a set of time instances, 
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TABLE I 

Parameters used for the simulations of the forced linear oscillator with stop. 





a 


c 


^max 


«(*) 


xo 


Xo 


3?max 




Example 1 


0.05 


0.9 


40tt 


20 cos(2.5t) 


11.36263 


31.40358 


14 


2.5 


Example 2 


0.95 


0.5 


4ir 


cos(t) 


-0.8 





-0.8 


1 



{*Jfe}fceN> computes a set of positions, z PS : {tfe} fcgN ->■ K, by: 

^ps(to) = £o, zps(*i) = %o + i h + y (u(0) - 2ax - w 2 x ), 
^ps(4+i) = -czps(tfc-i) + min{?/ps(tfc) ? (1 + c)x max }, 
l/ps(4) = ^-^(/i 2 M(t fe ) + (2-/iV)z PS (t fc )- ((l-c)-(l + c)a/i)zps(*fc-i)). 



(28) 



where t^+i = tk + h for each fceN. 

We illustrate the performance of our approach by considering the two examples described in Table [I] 
whose solutions, which are defined for all t G [0, t max ], can be computed analytically. The position 
component of the analytical trajectory of each example is plotted in Figs. |10b| and |10c[ The evaluation 
of t he pe rformance of our algorithm as described in Fig. [9] using p £ 



22, is shown in 



as in Definition 

Fig. |10d| To make our approach comparable to the PS Method, for A h we use a Runge-Kutta of order two 
which is called the midpoint method. We cannot use p £ to compare our discrete approximation algorithm 
to the PS method since the PS method does not compute the velocities of the hybrid system. Hence, we 
use the evaluation metric proposed in [|29ll which compares a numerically simulated position trajectory, 



^pos 



{tk} 



keN 



points {t k } keN D [0,/; max ] as follows: 



to the analytically computed position trajectory, x ana i y ti c : [(Mn 



— > R, at the sample 



/K^pos; -^analytic) max| |^pos(^fc) -^analytic (pk) I I {^fclfcgN ^ [^j^max]}- 



(29) 



The result of this comparison is illustrated in Fig.|10e Finally, the computation time on a 32 GB, 3.1 GHz 



Xeon processor computer for each of the examples as a function of the step-size and relaxation parameter 



is shown in Fig. lOf Notice in particular that we are able to achieve higher accuracy with respect to the p 
evaluation metric at much faster speeds. In Example 1, for step-sizes h < 10 _1 , our numerical simulation 
method is consistently more accurate by several orders of magnitude and generally several orders of 
magnitude faster than the PS method. In Example 2, using a step-size of approximately h = 10~ 2 and 
relaxation parameter e — 2 ■ 1CT 7 , our numerical simulation achieves a p value of approximately 1CT 4 
while taking approximately 0.1 seconds, whereas the PS method requires a step-size of h — 5 • 1CT 4 
which takes approximately 5 seconds in order to achieve the same level of accuracy. 



B. Navigation Benchmark for Hybrid System Verification 

Next, we illustrate the utility of our discrete approximation algorithm in Fig. [9] by considering three 
instances of a navigation benchmark proposed in If3~0il for hybrid system verification tools. The benchmark 
considers an object moving in the plane while following a set of desired velocities, Vd = (sin^) , cos^)), 
for j G {0, . . . , 7} where j is attributed to unit-sized squares according to a labeling map. Special symbols 
denoted "Goal" and "Obstacle" are reserved for a set of target and forbidden states, respectively. The 



labeling map for the three instances considered within this subsection are illustrated in Fig. 1 1 , where the 
label j in each cell refers to the desired velocity, target, or forbidden states. If the trajectory leaves the 
grid, the desired velocity is the velocity of the closest cell. 

The dynamics of the four dimensional state, (x,v) G IR 4 , are giv en by x (t) — v(t), and v(t) 



A(v(t)-v dj ), where A = ( 



1.2 
0.1 



0.1 
1.2 



) for the instances illustrated in Figs. 



11a 



and 



lib 



and A = [ _ 



-0.8 
-0.1 



-0.2 
-0.8 



) 



for the instance illustrated in Fig. |1 lc[ For each instance, we attempt to verify that for all trajectories 



beginning from a set of initial conditions there exists some finite time at which the "Goal" set is reached 
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(c) Infinite Instance 



Fig. 11. Navigation benchmark instances considered. Each instance map (symbols within each square) is shown with its desired velocity 
(drawn in each square with an arrow). Sample trajectories (drawn as lines beginning at filled in circles and ending at crosses) begin at some 
of the initial conditions which we attempt to verify. 



while avoiding the "Obstacle" set. We perform this verification by discretizing over the given set of initial 
conditions. 

For the instance illustrated in Fig. 11a, we select a set of initial conditions x G [0,1] x [0,1] and 
v G [0.1, 0.5] x [0.05, 0.25]. By choosing 10, 000 uniformly spaced points over the set of initial conditions, 
a step-size of 10 _1 , and relaxation size of 10~ 3 , we are able to verify this system in approximately 
100 seconds. For the instance illustrated in Fig. 1 lb we select a set of initial conditions x G [3, 4] x [3, 4] 
and v G [—1, 1] x [—1, 1]. This instance fails the verification task as trajectories are unable to reach the 
"Goal" set, but do not ever intersect with the "Obstacle" set. By choosing 10, 000 uniformly spaced points 
over the set of initial conditions, a step-size of 10 _1 , and relaxation size of 10~ 3 , we discover that for 



this system the task is not verifiable in approximately 85 seconds. For the instance illustrated in Fig. 11a 



we select a set of initial conditions x G [3,3.5] x [3,3.5] and v G {0.5} x [—0.5,0.5]. In this instance, 
verification is possible, but trajectories get close to the "Obstacle" set. By choosing 10, 000 uniformly 
spaced points over the set of initial conditions, a step-size of 10 _1 , and relaxation size of 10~ 3 , we are 
able to verify this system in approximately 210 seconds. 



C. Simultaneous Transitions in Models of Legged Locomotion 

As a terrestrial agent traverses an environment, its appendages intermittently contact the terrain. Since the 
equations governing the agent's motion change with each limb contact, the dynamics are naturally modeled 
by a controlled hybrid system with discrete modes corresponding to distinct contact configurations. 
Because the dynamics of dexterous manipulation are equivalent to that of legged locomotion [3TL such 
controlled hybrid systems model a broad and important class of dynamic interactions between an agent 
and environment. 

Legged animals commonly utilize gaits that, on average, involve the simultaneous transition of multiple 
limbs from aerial motion to ground contact 11321 . 11331 . Similarly, many multi-legged robots enforce simulta- 
neous leg touchdown via virtual constraints implemented algorithmically OH, 113511 or physical constraints 
implemented kinematically 11361 . [T3~71 . Trajectories modeling such gaits pass through the intersection of 
multiple transition surfaces in the corresponding controlled hybrid system models. Therefore simulation 
of this frequently-observed behavior requires a numerical integration scheme that can accommodate 
overlapping guards. Algorithm 11 has this capability, and to the best of our knowledge is the only existing 
algorithm possessing this property. We demonstrate this advanced capability using a pronking gait in a 
saggital-plane locomotion model. 

Fig. 12a| contains an extension of the "Passive RHex-runner" in [38] that allows pitching motion. A 
rigid body with mass m and moment-of-inertia I moves in the saggital plane under the influence of 
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(a) Schematic for the saggital-plane locomotion model 
with three mechanical degrees of freedom. 




(b) Projection of guards in (6, z) coordinates for transition from aerial 
domain D a to ground domain D g with parameters d — I — 1, 

ip — it/5. 



Fig. 12. Schematic and discrete mode diagram for the saggital-plane locomotion model. 
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Fig. 13. Snapshots of pronk at discrete transition times from initial condition [xo, zo, 0o, xo, &o, So) = (0,1.1,0,3.4,0,0), parameters 
(m, I , k, £, d, g,Tp) — (1, 1, 30, 1, 1, 9.81, vr/5), step size h = 10 -3 , relaxation parameter e = 10~ 2 (left). Same as before, but with 
6 = -0.4 (right). 



gravity g. Linear leg-springs are attached to the body via a frictionless pin joint located symmetrically at 
distance d/2 from the center-of-mass. The leg-springs are massless with linear stiffness k, rest length £, 
and make an angle ip with respect to the body while in the air. When a foot touches the ground it attaches 
via a frictionless pin joint, and it detaches when the leg extends to its rest length. 

A pronk is a gait wherein all legs touch down and lift off from the ground at the same time ll32l . 11331 . 
Due to symmetries in our model, motion with pitch angle 9 = for all time is invariant. Therefore periodic 
orbits for the spring-loaded inverted pendulum model in lf39ll correspond exactly to pronking gaits for 
our model. Fig. 12b contains a projection of the guards Gun, Gi a , r ), Gn g \, Gr r>g -\ in (9, z) coordinates for 
the transition from the aerial domain D a to the ground domain D g through left stance Di and right stance 
D r . The pronking trajectory is illustrated by a downward-pointing vertical arrow, and a nearby trajectory 
initialized with negative rotational velocity is illustrated by a dashed line. Fig. [13] contains snapshots from 
these simulations. 



The 6> = trajectory in Fig. 12b clearly demonstrates the need for a simulation algorithm that allows 
the intersection of multiple transition surfaces. We emphasize that our state-space metric was necessary 
to derive a convergent numerical approximation algorithm that accommodates this phenomena: since the 
discrete mode sequence differs for any pair of trajectories arbitrarily close to the 9 = execution that 
pass through the interior of Di and D r , respectively, a naive application of the trajectory-space metric 
in lfT5ll would yield a distance larger than unity between the pair. Consequently no numerical simulation 
algorithm can be shown to converge to the 9 = execution using a trajectory-space metric. As a final 
note to practitioners, we remark that our algorithm does not require a specialized mechanism to handle 
overlapping guards or control inputs: a single code will accurately simulate any orbitally stable execution 
of the hybrid system under investigation, dramatically simplifying practical implementation. 
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VI. Conclusion 

We developed an algorithm for the numerical simulation of controlled hybrid systems and proved the 
uniform convergence of our approximations to executions using a novel metrization of the controlled 
hybrid system's state space. The metric and the algorithm impose minimal assumptions on the hybrid 
system beyond those required to guarantee deterministic executions. Beyond their immediate utility, it 
is our conviction that these tools provide a foundation for formal analysis and computational controller 
synthesis in a broad class of CPS. 

There are many areas within the CPS community where we expect that our metrization and simulation 
framework will find fruitful application. For example, Girard and Pappas ll40l developed a family of 
approximate bisimulation metrics enabling comparison of entire CPS once a trajectory metric is provided 
for each particular system. We developed a general method to construct metrics on the state space and 
hence the space of trajectories for controlled hybrid systems, significantly extending the class of CPS that 
can be studied in this paradigm. Further, simulation provides a foundation for numerical tools including 
reachability-based controller synthesis iflTTl and numerical optimal control ll42ll . Current approaches to 
these problems require a fixed discrete mode sequence, yielding a computational complexity combinatorial 
in the number of discrete modes. We conjecture that our relaxed state-space metric enables generalizations 
of these algorithm which avoid combinatorial search by working in our continuous metric space. 

References 

[1] A. Nerode and W. Kohn, "Models for hybrid systems: Automata, topologies, controllability, observability," in Proceedings of the 

Workshop on the Theory of Hybrid Systems, 1993, pp. 317-356. 
[2] S. N. Simic, K. H. Johansson, J. Lygeros, and S. S. Sastry, "Towards a geometric theory of hybrid systems," Dynamics of Continuous 

Discrete and Impulsive Systems Series B: Applications & Algorithms, vol. 12, no. 5/6, pp. 649-688, 2005. 
[3] A. D. Ames and S. S. Sastry, "A homology theory for hybrid systems: Hybrid homology," in Proceedings of the 8th International 

Workshop on Hybrid Systems: Computation and Control, 2005, pp. 86-102. 
[4] K. H. Johansson, M. Egerstedt, J. Lygeros, and S. S. Sastry, "On the regularization of Zeno hybrid automata," Systems & Control 

Letters, vol. 38, no. 3, pp. 141-150, 1999. 
[5] L. Tavernini, "Differential automata and their discrete simulators," Nonlinear Analysis: Theory, Methods & Applications, vol. 11, no. 6, 

pp. 665-683, 1987. 

[6] D. Gokhman, "Topologies for hybrid solutions," Nonlinear Analysis: Hybrid Systems, vol. 2, no. 2, pp. 468-473, 2008. 
[7] D. Pollard, Convergence of Stochastic Processes, ser. Springer Series in Statistics. Springer, 1984. 

[8] C. Cai, R. Goebel, and A. R. Teel, "Relaxation results for hybrid inclusions," Set-Valued Analysis, vol. 16, pp. 733-757, 2008. 
[9] R. G. Sanfelice and A. R. Teel, "Dynamical properties of hybrid systems simulators," Automatica, vol. 46, no. 2, pp. 239-248, 2010. 
[10] M. Carver, "Efficient integration over discontinuities in ordinary differential equation simulations," Mathematics and Computers in 

Simulation, vol. 20, no. 3, pp. 190-196, 1978. 
[11] L. F. Shampine, I. Gladwell, and R. W. Brankin, "Reliable solution of special problems for ODEs," ACM Transactions on Mathematical 

Software, vol. 17, no. 1, pp. 11-25, 1991. 
[12] J. Guckenheimer and A. Nerode, "Simulation for hybrid systems and nonlinear control," in Proceedings of the 31st IEEE Conference 

on Decision and Control, 1992, pp. 2980-2981. 
[13] J. Esposito, V. Kumar, and G. J. Pappas, "Accurate event detection for simulating hybrid systems," in Proceedings of the 4th International 

Workshop on Hybrid Systems: Computation and Control, 2001, pp. 204—217. 
[14] L. Paoli and M. Schatzman, "A numerical scheme for impact problems II: The multidimensional case," SIAM journal on numerical 

analysis, vol. 40, no. 2, pp. 734-768, 2003. 
[15] L. Tavernini, "Generic asymptotic error estimates for the numerical simulation of hybrid systems," Nonlinear Analysis: Hybrid Systems, 

vol. 3, no. 2, pp. 108-123, 2009. 

[16] S. Burden, H. Gonzalez, R. Vasudevan, R. Bajcsy, and S. S. Sastry, "Numerical integration of hybrid dynamical systems via domain 

relaxation," in Proceedings of the 50th IEEE Conference on Decision and Control, 2011, pp. 3958-3965. 
[17] J. L. Kelley, General Topology, ser. Graduate Texts in Mathematics. Springer, 1955. 

[18] D. Burago, Y. Burago, and S. Ivanov, A Course in Metric Geometry, ser. Graduate Studies in Mathematics. American Mathematical 
Society, 2001. 

[19] M. Vidyasagar, Nonlinear Systems Analysis, 2nd ed., ser. Classics in Applied Mathematics. SIAM, 2002. 
[20] J. M. Lee, Introduction to Smooth Manifolds, ser. Graduate Texts in Mathematics. Springer, 2003. 

[21] V. I. Utkin, "Variable structure systems with sliding modes," IEEE Transactions on Automatic Control, vol. 22, no. 2, pp. 212-222, 
1977. 

[22] A. F. Filippov, Differential Equations with Discontinuous Righthand Sides. Kluwer Academic Publishers, 1988. 
[23] J. Zhang, K. H. Johansson, J. Lygeros, and S. S. Sastry, "Zeno hybrid systems," International Journal of Robust and Nonlinear Control, 
vol. 11, no. 5, pp. 435-451, 2001. 

[24] J. Lygeros, K. H. Johansson, S. N. Simic, J. Zhang, and S. S. Sastry, "Dynamical properties of hybrid automata," IEEE Transactions 
on Automatic Control, vol. 48, no. 1, pp. 2-17, 2003. 



21 



[25] E. Polak, Optimization: Algorithms and Consistent Approximation, ser. Applied Mathematical Sciences. Springer, 1997. 
[26] W. P. Ziemer, Weakly Differentiate Functions, ser. Graduate Texts in Mathematics. Springer, 1989. 

[27] R. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems, 

ser. Classics in Applied Mathematics. Society for Industrial and Applied Mathematics, 2007. 
[28] W. Rudin, Principles of Mathematical Analysis. McGraw-Hill, 1964. 

[29] O. Janin and C. Lamarque, "Comparison of several numerical methods for mechanical systems with impacts," International Journal 

for Numerical Methods in Engineering, vol. 51, no. 9, pp. 1101-1132, 2001. 
[30] A. Fehnker and F. Ivancic, "Benchmarks for hybrid systems verification," in Proceedings of the 7th International Workshop on Hybrid 

Systems: Computation and Control, 2004, pp. 326-341. 
[31] A. M. Johnson, G. C. Haynes, and D. E. Koditschek, "Standing self-manipulation for a legged robot," in Proceedings of the IEEE 

International Conference on Intelligent Robots and Systems, 2012. 
[32] R. M. Alexander, "The gaits of bipedal and quadrupedal animals," The International Journal of Robotics Research, vol. 3, no. 2, pp. 

49-59, 1984. 

[33] M. Golubitsky, I. Stewart, P. L. Buono, and J. J. Collins, "Symmetry in locomotor central pattern generators and animal gaits," Nature, 

vol. 401, no. 6754, pp. 693-695, 1999. 
[34] M. Raibert, M. Chepponis, and H. Brown Jr., "Running on four legs as though they were one," IEEE Journal of Robotics and Automation, 

vol. 2, no. 2, pp. 70-82, 1986. 

[35] U. Saranli, M. Buehler, and D. E. Koditschek, "RHex: A simple and highly mobile hexapod robot," The International Journal of 

Robotics Research, vol. 20, no. 7, pp. 616-631, 2001. 
[36] S. Kim, J. Clark, and M. Cutkosky, "iSprawl: Design and tuning for high-speed autonomous open-loop running," The International 

Journal of Robotics Research, vol. 25, no. 9, pp. 903-912, 2006. 
[37] A. Hoover, S. Burden, X. Fu, S. S. Sastry, and R. Fearing, "Bio-inspired design and dynamic maneuverability of a minimally actuated 

six-legged robot," in Proceedings of the 3rd IEEE International Conference on Biomedical Robotics and Biomechatronics, 2010, pp. 

869-876. 

[38] J. Seipel and P. Holmes, "Three-dimensional translational dynamics and stability of multi-legged runners," The International Journal 

of Robotics Research, vol. 25, no. 9, pp. 889-902, 2006. 
[39] R. M. Ghigliazza, R. Altendorfer, P. Holmes, and D. E. Koditschek, "A simply stabilized running model," SIAM Journal on Applied 

Dynamical Systems, vol. 2, no. 2, pp. 187-218, 2003. 
[40] A. Girard and G. J. Pappas, "Approximation metrics for discrete and continuous systems," IEEE Transactions on Automatic Control, 

vol. 52, no. 5, pp. 782-798, 2007. 

[41] J. Ding, J. H. Gillula, H. Huang, M. P. Vitus, W. Zhang, and C. J. Tomlin, "Hybrid systems in robotics," IEEE Robotics & Automation 

Magazine, vol. 18, no. 3, pp. 33-43, 2011. 
[42] E. Westervelt, J. Grizzle, and D. E. Koditschek, "Hybrid zero dynamics of planar biped walkers," IEEE Transactions on Automatic 

Control, vol. 48, no. 1, pp. 42-56, 2003. 



